Constant-temperature molecular-dynamics algorithms for mixed hard-core/continuous 
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We present a set of second-order, time-reversible algorithms for the isothermal (NVT) molecular- 
dynamics (MD) simulation of systems with mixed hard-core/continuous potentials. The methods are 
generated by combining real-time Nose thermostats with our previously developed Collision Verlet 
algorithm [Mol. Phys. 98, 309 (1999)] for constant energy MD simulation of such systems. In all 
we present 5 methods, one based on the Nose-Hoover [Phys. Rev. A 31, 1695 (1985)] equations of 
motion and four based on the Nose-Poincare [J.Comp.Phys., 151 114 (1999)] real-time formulation 
of Nose dynamics. The methods are tested using a system of hard spheres with attractive tails and 
all correctly reproduce a canonical distribution of instantaneous temperature. The Nose-Hoover 
based method and two of the Nose-Poincare methods are shown to have good energy conservation 
in long simulations. 
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I. INTRODUCTION 

Algorithms for molecular-dynamics simulation can be 
generally be divided into two distinct classes depending 
upon the nature of the potential For systems with 
continuously differentiable potentials, the trajectory is 
generated through the numerical integration of the equa- 
tions of motion - a coupled set of differential equations - 
typically with a fixed time step. At the other end of the 
spectrum are methods for discontinuous potentials, such 
as hard spheres or the square-well potential. Such algo- 
rithms are event driven in that the system is advanced 
ballistically between " collisions" , which are then resolved 
exactly. There exist, however, model interaction poten- 
tials of theoretical and practical importance that are hy- 
brids of continuous and discontinuous potentials - for ex- 
ample, the restricted primitive model for electrolyte solu- 
tions or the dipolar hard-sphere model of polar fluids. To 
date, simulation studies for such systems have primarily 
been restricted to Monte Carlo studies due to the lack of 
a viable molecular-dynamics (MD) algorithm. To rem- 
edy this, we have recently introduced a new molecular- 
dynamics method for such systems The algorithm, 
referred to as Collision Verlet, has good energy conserva- 
tion and is far more stable over long time simulation than 
previous integrators for hybrid continuous/discontinous 
systems. The Collision Verlet algorithm was formulated 
as a constant energy simulation method, which generates 
configurations from a microcanonical (NVE) distribu- 
tion. However, to mimic experimental conditions most 
modern simulations are run under isothermal (NVT) or 
isothermal/isobaric (NPT) conditions. In this work, we 
introduce and evaluate several reformulations of Collision 
Verlet to generate trajectories whose phase space points 
are canonically (isothermally) distributed. 



The NVT (isothermal) Collision Verlet algorithms de- 
veloped here are all based on the extended Hamilitonian 
of Nose|^, which is a standard technique for generating 
canonical trajectories for the simulation of systems with 
continuous interaction potentials. In the Nose approach, 
the phase space of the system is augmented by the in- 
troduction of an auxilliary variable s and its conjugate 
momentum tt (with "mass" Q). For a system with a 
potential V, the Nose extended Hamiltonian is 
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where pi is the momentum conjugate to the position qi 
and is related to the actual momentum , pi, by the rela- 
tion Pi = Pi/ s, and the parameter g — Nf + 1, where Nf 
is the number of degrees of freedom of the system. With 
this choice oi g, it can be readily shown assuming er- 
godicity, that constant energy (microcanonical) dynamics 
generated by the Nose Hamiltonian produces a canonical 
(constant temperature) distribution in the reduced phase 
space {p/s,q}. 

The generation of phase space configurations dis- 
tributed in the canonical ensemble within the Nose dy- 
namical scheme is accomplished by a dynamical rescaling 
of time. The real time of the simulation, t, is related to 
the Nose time, r, by the transformation ^ = s. Since 
numerical integration methods generally operate with a 
fixed time step, the transformation to real time generates 
a nonuniform grid of time points Q, which is inconvenient 
for the calculation of system averages. To remedy this, 
two schemes have been developed to produce equations 
of motion for Nose dynamics that generate trajectories 
directly in real time. By applying time and coordinate 
transformations directly to the Nose equations of mo- 
tion Hoover derived a set of real-time equations of 
motion for Nose dynamics, defining the so-called Nose- 
Hoover method. This approach has become the most 
widely isothermal simulation method, but has a draw- 
back in that the coordinate transformation used is not 
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canonical and the Nose-Hoover equations of motion are 
non-Hamiltonian in structure, precluding the use of sym- 
plectic integration schemes 0. In an alternate approach 
Bond, Leimkuhler and LairdQ apply a Poincare time 
transformation to the Nose Hamiltonian to give the so- 
called Nose-Poincare Hamiltonian, from which real-time, 
fully Hamiltonian equations of motion for Nose dynamics 
are generated. 

In this work we present constant temperature simu- 
lation methods for mixed continuous/discontinuous in- 
teraction potentials generated by adapting the Collision 
Verlet method within both the Nose-Hoover and Nose- 
Poincare schemes. In the next section we briefly review 
the standard Collision Verlet algorithm[|j followed by the 
introduction of the Nose-Hoover Collision Verlet (NHCV) 
and Nose-Poincare Collision Verlet (NPCV) algorithms 
in Sections 3 and 4, respectively. The algorithms are 
evaluated in Section 5 through numerical experiments on 
a model system. In section 6, we conclude. 

II. THE COLLISION VERLET ALGORITHM 

In this section we review the Collision Verlet Q al- 
gorithm for the numerical integration of the dynamics 
of systems with mixed continuous/discontinuous inter- 
action potentials. We consider N particles interacting 
through a continous potential plus a hard core, assumed 
here to be spherical. To facilitate the construction of nu- 
merical methods, it is useful to describe the dynamics of 
the system within a Hamiltonian format, but for a sys- 
tem with a discontinuous potential the construction of 
a Hamiltonian as the generator of the dynamical equa- 
tions of motion is problematic. In this work, we observe 
that the hard sphere interaction potential, 14s ({q}) typ- 
ically can be approximated to any degree of accuracy by 
a sequence of steeply repulsive continuous functions. In 
this sense, the energy function Ti. of the mixed system 
is refered to here as a pseudo-Hamiltonian. Here the 
pseudo-Hamiltonian is given by 

H = T{p) + VU{<l}) + Vc{{q}), (2) 

2 

where the kinetic energy T(p) = Vhsiin}) is 

the hard sphere potential and Vc({q}) is a continuously 
differentiable potential energy function, that we assume 
to be pairwise additive, that is, 

i j>i 

where Vc is a pair potiential, qij is the distance between 
two particles indexed by i and j, and the sum is over all 
pairs of particles. 

The Collision Verlet algorithm is based on the splitting 
of the continuous pair potential, Vc{q), into a short range 
part, vi{q), and a long range part, V2{q), according to 

Vc{q) = vi{q) + V2{q) (3) 



The potential splitting is rendered so that the force due 
to the long-range part of the potential vanishes at the 
hard-sphere contact distance (i.e. W2(cr) — 0). This 
form of the potential splitting is necessary for the con- 
struction of a second-order method - For the motivation 
and specific details of this splitting technique the reader 
is referred to reference Q. The pseudo-Hamiltonian is 
then split accordingly. For generality, let consider Ti. 
as a pseudo-Hamiltonian of any given mixed impulsive- 
continuous system. Next, we partition Ti, in the following 
way: 

W = 7^l + H2 , (4) 

where Tii includes the kinetic energy, the hard sphere 
potential, Vhs, and the short range potential, Vi; Ti2 
must include the long range potential, V2. A Trotter 
factorization then gives the following approximation 
for the dynamical flow map, (/)>^(t), defined as the oper- 
ator (associated with the Hamiltonian Ti) that advances 
the phase space configuration a time r into the future, 

T T 

Mt) = <I)-H2 {^)4>'Hi {r)(f>'H2 ( 2 ) (5) 

Since only contains the long-range potential, the flow 
map 0-^2 can be constructed exactly. The flow map cor- 
responding to Til is approximated in the following way 

1=1 

where Uc is the number of hard-sphere collisions during 

(c) 

the time step h, is the time between each collision 
(with t['^'^ being measured from the beginning of the time 

(c) 

step until the first collision and r„ Yi measured from the 
last collision to the end of the time step so that J^^li^ = 
r), and (t)Vhs is flow map for an instantaneous hard- 
sphere collision. Finally, the flow map for the motion of 
the particle between collisions is approximated using the 
Stomer- Verlet algorithm generated by a further Trotter 
factorization 

<Pt+vAt) ~ <l>vA^)<l>T{T)(j)vA^) ■ (7) 

The most CPU intensive part of the Collision Verlet 
algorithm is the determination of the time to next colli- 
sion Tc- The collision condition for two particles i and j 
can be written as 

lk.(rc)-<Z,(r.)f -a^-O. (8) 

Since the trajectories between collisions are approxi- 
mated within the Collision Verlet scheme by quadratic 
equations, the collision condition (^ is a quartic equa- 
tion. To ensure that all collisions are resolved correctly, 
it is necessary to accurately resolve the smallest positive 
root to this quartic equation. This is not a trivial prob- 
lem as the root becomes increasingly unstable as smaller 
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time steps are used (i.e., when the time to cohision is 
small). To increase efficiency and accuracy of the com- 
putation, we employed in all the simulations in this paper 
a root finding method based on Cauchy indices The 
details of the collision-time calculation are given in the 
Appendix. 



III. COLLISION VERLET WITH A 
NOSE-HOOVER THERMOSTAT 

As discussed in the introduction, the Nose-Hoover 
method for isothermal molecular-dynamics simulation is 
generated by applying time and coordinate transforma- 
tions to the equations of motion generated by the Nose 
Hamiltonian (Eq. |^), which are 
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dr dqi 



Conversion to real time, t, is accomplished through the 
following transformations 



p dr 

s at 



(11) 



In addition. Hoover simplified the resulting equations of 
motion by introducing a further variable tranformation 



?7 = In s ^ = ?? 



(12) 



resulting in the so-called Nose-Hoover equations of mo- 
tion: 

* = — , Pi = -^V{q) -pi^, (13) 



with the earlier literature, we write the flow map in terms 
of a Liouville operator, £, as follows 



(r) = 



(16) 



The Liouville operator corresponding to the Nose-Hoover 
equations of motion above is 



dpi ^ dqi dp, 

L 

,2 \ Q 



(17) 



where we have explicitly included a hard-sphere term. 

To get a reversible method for the Nose- Hoover method 
with mixed potentials, the above Liouville operator is 
split in the following way: 



with 



£1 = Lhs + 



Pi d 

nil dqi ^ dqi ' dp. 



dqi dpi 



(18) 



(19) 
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and 
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A Trotter factorization is now applied to this splitting. 



(22) 
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These equations of motion can be shown to generate 
configurations distributed according to an isothermal 
(canonical) distribution as long as the system is ergodic 
and g — Nf, the number of degrees of freedom. 
Since the coordinate transformation is non-canonical, the 
equations of motion are not derivable from a Hamilto- 
nian, however a conserved energy does exist and is given 
by 



^ 2m, 
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In order to simplify the construction of splitting meth- 
ods for this non-Hamiltonian system and to make contact 



The operator e'~^'^ is approximated using the Collision 
Verlet method described in the previous section - see 
Eq. ^. The solution of the operator e'^^'^/^ is straight- 
forward. To find the solution of the operator e'''^'^/^,i.e, 
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we further split £3. That is, 



with 
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The corresponding Trotter factorization of this spUtting 
is 



(27) 



The solution of the operator e 3 '^/^ is straightforward. 

^(1) 

The operator e 3 is solve from a further splitting. The 
solution of the operator e'-^'^l'^ gives 




- "Tin + 2^)1+1/2, 



Pi,n+1 = Pi,' 



1 - ^en+l/2/4 

l + re„+i/2/4' 



(29) 



(30) 



6..r^W + ^(E^-.^^)- (31) 

The algorithm is tested in Section 5 for a system of 
hard spheres with inverse-sixth-power attractive tails. 

Certainly, the Liouville operator splitting used above 
is not the only possible method. For example, another 
splitting is 



£ = £1 + £2 



(32) 



with 



^1 = E -7^ + ^'^^ - E 7^^i('^)7^' (33) 

^ mi dqi ^ aqi dpi 



and 

£2 = ~y^P^i-^ - ^^2(9)^ 
^ dpi dqi dpi 

can be used. Using a Trotter factorization gives 



(34) 



(35) 



for systems with continuous potentials, can be shown 
to enhance long-term stability]^. Recently, Bond, 
Leimkuhler, and Laird]^ have proposed a new real-time, 
but fully Hamiltonian, formulation of the Nose constant- 
temperature dynamics. This is accomplished by perform- 
ing a time transformation, not to the Nose equations of 
motion as with Nose-Hoover, but directly to the Hamil- 
tonian using a Poincare time transformation, as follows: 



y-NP — siTi-Nose — Wo), 



(36) 



where Hq is the initial value of Ti. Nose ■ Combining equa- 
tions (|l]) and (|3^) the Nose-Poincare thermostat Hamil- 
tonian of a physical system consisting of N particles is 
expressed as following 



Hnp = s 



(E^ + ^^('^) + g+5^^i--^o). 

(37) 

In order to sample the correct canonical distribution, 
the constant g is taken to be the number of degrees of 
freedom|7|, g = Nf. The equations of motion are 



Pi . TT 

miS Q 



(38) 



dqi rn 



miS'^ 



gkT - AH, (39) 



Pi 



2m, s' 



2 ■ VM + ^+gkTlns-Ho. (40) 



Note that, the exact solution to Nose-Poincare equations 
of motion generates trajectories that are identical to that 
generated by the Nose-Hoover scheme, exactly solved. It 
is in the construction of approximate numerical methods 
that these two approaches differ. 

For the present case, we write the Nose-Poincare ther- 
mostat pseudo-Hamiltonian (see Sect. 2) for a mixed 
hard-core/continuous potentials system 



n 



NP 



-gkT In s — Ha) 



(41) 



There are a variety of ways in which one can construct 
numerical integration algorithms using this Hamiltonian. 
To this end, we first consider two ways of splitting the 
overal NP Hamiltonian:: 

Splitting I 



IV. COLLISION VERLET WITH A 
NOSE-POINCARE THERMOSTAT 

The Nose- Hoover formulation of constant-temperature 
dynamics is non-Hamiltonian in structure, thereby pre- 
venting the use of symplectic integration schemes, which. 



^i = ME^ + ^''^('?) + ^i('?) 



H2 = s[ 1/2(9) 



gkT\ns-Ho) (42) 



2Q 



(43) 
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Splitting II 

^1 = ^fE^ + ^''^(9) + ^i(<z)-Woj (44) 

(45) 



^2 = s\^V2iq) + — +gkTlns 



A Trotter factorization of the flow map (Eq. ^ is applied 
to each splitting. To approximate the flow map gener- 
ated by Til , we employ the Collision Verlet Scheme given 
in Eq. || to integrate the system from collision to collision 
under the influence of the short-range potential. Since s 
is a constant in the dynamics generated by Tii in both 
splittings, the Stormer- Verlet algorithm can be used to 
integrate the trajectory between collisions, with the colli- 
sion time being calculated as described in the Appendix. 
For splitting I, Stormer- Verlet gives 



Pi,n+l/2 



Pi,n+l/2 - T^^n^Viyqn 

2 dm 



7r„+i/2 + 2 



1 ( Pi,n+l/2 
n; Pi,n+l/2i 



Qi,n+1 



_Pi,n+l/2 



TTn+l = 7r„+i/2 + 



T 



^ — ^ m : 



Pi,n+l/2 



-AH (g„+i,Pi,„+i/2, s„) 

Pi,n+1 = Pi,n + l/2 ' 



^Sn^V"i((7„+i) 

2 dqi 



(46) 

(47) 
(48) 

(49) 
(50) 



The equations for Splitting II can be similarly generated. 

In both Splittings I and II the integration of H2 is com- 
plicated by the presence of both s and its conjugate mo- 
mentum TT, but here we consider two possible approaches: 

7^2 Integration Method 1: Since the Hamilto- 
nian here is non-separable, the Generalized 
Leapfrog 0, |l^ scheme, a fully symplectic 
extension of the Stormer- Verlet algorithm for 
non-seperable Hamiltonians, can be used. The 
integration for Splitting I for timestep r is 



Pi,n+l/2 — Pi,n 



^S„-^V2((?rO 



(51) 



7r„+i/2 = T^sm - ^ {gkT + AH2 (g™, s„, 7r„+i/2)) (52) 



, T N T^n+1/2 
Sn+1 = Sn + - [Sn + S„+i) — , 



(53) 



'ri+l 



^«+i/2 - 2 (■9^'^ + Ai?2 (g™, S„+l, 7r„+i/2)) 



Pi,n+1 = Pi. 



-1/2 - ■7;^n+l^V2{qn) 

2 dqi 



(54) 



(55) 



The above integration is explicit. Eq. ^ requires 
the solution of a scalar quadratic equation for 
7r„_|_i/2- Details of how to solve this equation with- 
out involving subtractive cancellation can be found 
in Ref . |0] . The application of Method 1 for the I-L2 
in Splitting II is similar and straightforward. 

7^2 Integration Method 2: Instead of using General- 
ized Leapfrog, we employ a splitting of 7^2 



(1) 



For Splitting I, we use 



H 



(2) 



^(1) _ STT 



n 



(2) 



sV2{q) 



(56) 

(57) 
(58) 



Since no conjugate pair appears in 7^2^^ its dynam- 
ics for a timestep r is straightforward 

d 



Pi,n-\-l Pi.n 
T^n+l = 7r„ - TV2{qn) 



TSn^V2iqn) 

dqi 



(59) 
(60) 



Only equations involving variables p and tt are 
shown above because q and s are constants of mo- 
tion. 

The solution of the dynamics of 7^2^ involves a 
conjugate pair s and tt, but it can be solved ex- 
actly |ll]. Thus the time evolution of 7^2^-' for the 
timestep t is 



Sn+l 



1 + 

2Q 



1 



— r 

2Q 



(61) 



(62) 



Here, it is q, and p that are constants of motion. 
Again, the application of Method 2 for Splitting II 
is similar and straightforward. 

Combining the two overall splittings for the NP Hamil- 
tonian with the two methods for integrating 7^2, gives 
a total of 4 proposed algorithms for the Nose-Poincare 
Coflision- Verlet (NPCV) method. These are 

• NPCVl: Splitting I + 7^2 integration method 1 

• NPCV2: Splitting I +712 integration method 2 
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FIG. 1: order of accuracy of the NHCV algorithm and NPCV 
algorithms 1 to 4. Comparison is made with a line of order 2. 



• NPCV3: Splitting II +7^2 integration method 1 

• NPCV4: Splitting II +H2 integration method 2 

In the next section we test these four algorithms for a 
model system and compare them with each other and 
with the Nose-Hoover Collision Verlet (NHCV) method 
outlined in the previous section. 



V. NUMERICAL EXPERIMENTS ON A 
MODEL POTENTIAL 

We test the various algorithms for NVT Collision Ver- 
let proposed in this paper using a system of hard-spheres 
with an attractive inverse-sixth-power continuous poten- 
tial, 



.1 



(63) 



where a is the hard-sphere diameter. The potential is 
truncated at the distance Qc = 2.5a and, to ensure its 
continuity, it is shifted and smoothed so that potential 
and the force vanish beyond the cutoff distance. We split 
the above potential into short and long-range parts, as 
prescribed in Ref.[^, with qi and (72 as input parameters. 

The MD simulations were carried out on systems of 
N — 500 particles. A system of reduced units was cho- 
sen so that all quantities are dimensionless: as units of 
distance and energy we used the potential parameters 
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FIG. 2: Energy conservation in a long simulation run (10 
time steps) for NPCV algorithms 1 to 4. 
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FIG. 3: Energy versus time in a long simulation run (10^) 
using the NHCV and NPCVl algorithms 
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FIG. 4: Instantaneous temperature distributions for the 
NPCV algorithms 1 to 4. In each, the exact canonical dis- 
tribution is shown as a solid line. 
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FIG. 5: Instantaneous temperature distribution for the 
NHCV simulations (circles). The exact canonical distribu- 
tion is shown as a solid line. 



a and e, respectively, and the mass of one atom as the 
unit mass. The unit of time is (mcr^/e)"'^/^. An asterisk 
superscript indicates reduced units. In all simulations 
the density was p* = pa^ = 0.7 with reduced temper- 
ature T* ~ kT / e = 1.5. In addition, a cubic box with 
periodic boundary conditions was used. In improve effi- 
ciency, neig hbor (Verlet) lists § were used for the evalu- 
tion of the short range force, the long range force, and the 
collision times. In all of our simulations, we set g = Nf 
with Nf = 3 (TV — 1) to correct for the fact that in a 
molecular-dynamics simulation the total linear momen- 
tum is conserved p2|. Each run has was started form an 
initial configuration produced after an equilibration run 
of 200,000 time steps (with r* = 0.001) starting from an 
fee (face-centered-cube) lattice with the particle veloci- 
ties chosen from a Boltzmann distributuion at T* = 1.5. 
The initial values of the extended variables in all of the 
numerical experiments are set to be so = 1 and psfl = 
in the case of the Nose-Poincare thermostat methods. In 
the case of the Nose-Hoover method, the initial values of 
the extended variables are thus 770 = and ^0 = 0. 

In order to compare the short time accuracy of the 
methods and verify that each one exhibits second-order 
global error, we show in Figure 1 a log-log plot of the 
maximum energy error for a run of total length t* = 
12 for each method as a function of time step, t. For 
comparison, a line of slope 2 is plotted to show that the 
global error for each method is second order, as required. 
In these runs the thermostat mass Q was set to 1.0. Note 
that, due to the discontinuous nature of the dynamics, 
the second order global error is not simply a consequence 
of the time-reversibility of the algorithms, but it also a 
direct result of the particular potential splitting we have 
chosen|2|. From Figure 1 we see that for short runs, 
the Nose-Hoover based method has the smallest error 
constant. 

For molecular-dynamics simulation the stability during 
long runs is more important that the order or short-term 
behavior of the algorithm. To test these we plot the 
energy trajectory, 5E — E[t) — E{t = 0), versus time 
for each of our methods using 10^ time steps of length 
T* = 5 X 10~^ (total time 5000). Figure 2 shows this plot 
for each of the 4 Nose-Poincare based methods discussed 
in the previous section. For this system, NPCV methods 
2 and 3 exhibit significant drift whereas methods 1 and 
4 are more stable for long time trajectories. The same 
plot for the Nose-Hoover method presented in section 3 
is shown in Figure 3 with the plot for NPCV method 1 
shown for comparison. The NPCV method 1 has slightly 
better energy conservation for this system than the Nose- 
Hoover Collision Verlet algorithm, which is comparable 
to NPCV method 4, but the differences are small and 
could change depending on the system. 

The algorithms presented here are designed to give a 
canonical distribution of phase space points. A useful 
check of this is to examine the distribution of instan- 
taneous temperature (as defined for a system with zero 
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total momentum) 

f = = V ^ (64) 

3(iV - 1) ^ 2m ^ ' 

A canonical distribution in momenta requires that this 
quatitiy be Gaussian distributed about the target tem- 
perature T with a variance of -^^-ypj- In Figure 4 is 
plotted the temperature distributions for the 4 NPCV 
algorithms using a thermostat mass of 10 measured dur- 
ing runs of 270,000 time steps (r* = 5 x 10~^) after 
equilibration. Figure 5 shows the same quantity for the 
Nose-Hoover Collision Verlet method. Comparison with 
the theoretical distribution, shown as a solid line in each 
plot, indicates that the canonical distribution is well re- 
produced by all proposed algorithms. 



VI. CONCLUSION 

In this work we have developed several algorithms, 
based on the extended Hamiltonian thermostat of Nose, 
to perform constant temperature (NVT) molecular- 
dynamics simulations of systems with mixed hard- 
core/continuous potentials. The methods are extentions 
of our recently developed Collision Verlet method for 
constant energy (NVE) MD simulation of such systems. 
These new methods, to our knowledge, represent the first 
viable canonical molecular-dynamics simulation methods 
for hybrid discontinous/continuous potentials. 

Specifically, five new algorithms have been presented 
and tested. The first algorithm, the Nose-Hoover Colli- 
sion Verlet (NHCV) algorithm, is based on application 
of the Nose- Hoover thermostat to the Collision Verlet 
scheme. The other 4 algorithms presented are based on 
the Nose-Poincare formulation of real-time Nose dynam- 
ics. These Nose-Poincare Collision Verlet methods differ 
from one another in the details of the numerical scheme 
used to integrate the equations of motion. All meth- 
ods were shown to give second-order global error in test 
simulation with the NHCV method having the smallest 
error constant for short-time simulations. The NHCV 
algorithm and two of the presented NPCV algorithms 
(NPCVl and NPCV4) were found to exhibit good sta- 
bility in long time simulations involving 500 hard-sphere 
particles with attractive inverse-sixth-power tails. In ad- 
dition, all methods were shown to correctly reproduce 
the canonical distribution of instantaneous temperature 
(kinetic energy). Note that, if the continuous potential 
is set to zero, the presented methods also provide a way 
of performing canonical, as opposed to isokinetic, hard- 
sphere molecular-dynamics simulations. 
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APPENDIX A: CALCULATION OF TIME TO 
NEXT COLLISION 

In this appendix we address the issue of the collision 
time calculation for mixed hard-core/continuous poten- 
tials systems. The quartic equation for the collision con- 
dition (Eq. ^), is solved for all pairs of particles and the 
smallest positive root is located as the time to the next 
collision. For mixed hard-core/continuous potentials sys- 
tems, this is time-consuming operations since collision 
times for all pairs must be recalculated after each colli- 
sion. In addition, Eq. ^ is quartic and difficult to solve. 
As we said in section 2, the quartic equation must be 
solved accurately to give the nearest root to zero in or- 
der to make sure that no collisions are missed. 

In ref.Q, we employed Laguerre's method [p^ for colli- 
sion time calculation for mixed hard-core/continuous po- 
tentials systems. The method is sufficient for all but the 
very smallest timesteps studied. But the method turns 
out to be very slow. This because for any given time 
interval and pair of particles, all the four complex roots 
need to be calculated. Also Laguerre's method deals with 
complex arithmetic. In this appendix, we propose a time 
saving collision time calculation method for collision ver- 
let. This method is based on a Cauchy indices of a Sturm 
sequence of a real polynomial in a real interval. 

The Cauchy index is an integer that can be associated 
with any real rational function and any interval whose 
end points are not the function poles. Let r be a ratio- 
nal function. The Cauchy index, /fr(a;), of r for the 
interval [a, P] is by definition the number of jumps of the 
function r from -|-oo to —oo on the interval [a,/?]. The 
Cauchy index can be calculated for any real polynomial 
that forms a Sturm sequence, {fo, fi, fm}, for the 
interval [a,/?]. The definition of the Sturm sequence of 
a real polynomials can be found in ref.Q. The connec- 
tion between the Cauchy index and the number of sign 
changes, v{x) for arbitrary real x, in the numerical se- 
quence ,{/o, /i, /m}, is given by the following result 
due to Sturm 0. 



Theorem 1 Let the real polynomials ,{/o, /i, /m} 
form a Sturm sequence for the interval [a,P], a < (3. 



9 



Then 



7o 



v{a) — v{P). 



(Al) 



Using this theorem we can write the number of real 
roots for a given polynomial p in any real interval [a, (3] 
in terms of the Cauchy index 



Po 



(A2) 



of the sequence {pk}-, generated by the Euclidean 
algorithm Q using the starting polynomials po -^ViPi 
p', with p' being the first derivative of the polynomial p. 
The elements of the rest of the sequence are linked by 
the relations 



Pi{x) 



qi{x)pi{x) 
q2{x)p2{x) 



-P2{x), 
-Psix), 



(A3) 
(A4) 



use the Sturm-Habicht pseudodivisions subresultant 
(PRS) method Ip^. The members of the polynomial re- 
mainder sequence pi{x),p2{x),p3{x), ...,ph{x) 

lc[Pt+i{x)]"'~"''^'-'^^Pi{x) = pi+i{x)qi{x) - (3iPi+2{x), 

(A7) 



deg[pi+2ix)] < deg[pi+i{x) 



(A8) 



where i = l,2,..,/i — 1, for some h, rii = deg[pi{x)], 
and lc[pi{x)] is the leading coeficient of pi. The different 
values of (3i are 



f3i 

H2 



(A9) 

i = 2,3,...,/i- 1, (AlO) 
{iMxW'-'^^ (All) 

{/,[K(x)]}"'--"'i?,^~/"'-^-"'\ 

i = i,...,h-l (A12) 



Pk^i{x) = qk{x)pk{x) ~pk+i{x), 



Pr, 



v{x) = qm{x)pm{x). 



(A5) 



(A6) 



Let 



The Euclidean algorithm also furnishes information 
about the multiplicity of the zeros, xq is a zero of mul- 
tiplicity fc of p if and only if it is a zero of multiplicity 
fc — 1 of Pm- We are now able to develop a collision time 
calculation method for Collision Verlet. 

From the above, the first step for Collision Verlet colli- 
sion time calculation is to determine in a given time inter- 
val the number of real roots by calculating the Cauchy 
index for the time interval. This means that we need 
an algorithm for polynomial division. The main prob- 
lem with polynomials division is that the bitlenght of 
coefficients in the sequence can increase dramatically 
and also, because we are dividing, in some cases the 
denominator can vanish. To solve this problem, we 



P{x) 



ax 



bx-^ 



cx 



■ dx + 



(A13) 



be the quartic polynomial obtained from the collision 
condition of eq. (||), and {pi,P2,P3,P4:,P5} i ts S turm- 
Habitch sequence determined by using eq. A7. We 



now determine the number of real roots of the equation 
p{t) = in a given time interval by calculating its Cauchy 
index, Eq. A2. If there is only one root, then we use 
Newton-Raphson method (isj to approximate the root. 
If there is more than one root, then we combine bisec- 
tion method and root counting method to isolate 
the time interval containing the smalest root. 

This method for solving for the shortest collision time 
is quite efficient giving a factor of 20 speed-up from our 
previous simulations using the Laguerre method , pri- 
marily because we no longer calculate all four roots of 
the quadratic equation and avoid complex arithmetic. 
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